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We investigate numerically the dynamical behavior of the mean field 3-spin spin glass model: we 
study equilibrium dynamics, and compute equilibrium time scales as a function of the system size 
V . We find that for increasing volumes the time scales r increase like In r oc V . We also present an 
accurate study of the equilibrium static properties of the system. 

PACS numbers: PACS numbers: 75.50.Lk, 75.10.Nr, 75.40.Gb 
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INTRODUCTION 

The glassy state is a state of matter ubiquitous in na- 
ture. It enjoys a number of intriguing properties: among 
those is the amazingly slow dynamics, common to both 
structural glasses and spin glasses. It is surely not well 
understood, and during the last years it has been the 
subject of numerous investigations, both theoretical and 
experimental. 

One interesting approach, first used for studying spin 
glasses in is based on (optimized) Monte Carlo nu- 
merical simulations. One considers samples that can be 
brought to equilibrium (in order to do that the lattice size 
cannot be too large) and studies numerically the equilib- 
rium dynamics of the system (this is typically not what 
happens with real disordered systems of this type, that 
are forever out of equilibrium). 

The method was first used to investigate the dynamics 
of the mean field Sherrington-Kirkpatrick (SK) model 
it turns out 123 that all the relevant time scales r a 
of the model grow for diverging lattice sizes V according 
to the scaling law ln(r Q ) oc U 1 / 3 . A single time scale 
controls the critical behavior of the model: the time scale 
governing the mode leading from A, one of the many 
equilibrium states of the mean field theory, to a different 
state B (not related to A by a global spin flip) has the 
same scaling law that the time scale controlling the mode 
related to a Z 2 global inversion (leading from a state to 
its Zi symmetric). 

We will deal here with the case of the mean field, fully 
connected p-spin spin glass, with p = 3, and try to de- 
termine the relevant time scales. Recently Ioffe, Lopatin 
and Sherrington have developed a theoretical anal- 
ysis of the problem: their approach suggests that for 
this class of models the relevant time scale behaves as 
ln(r Q ) oc V (see also the work of 4] for a computa- 
tion of relaxation times in spin systems with disordered 
quenched couplings). Here we try to analyze the model 
to uncover the true asymptotic behavior by using the 



numerical approach of 0, 0] ■ 

The fully connected p-spin model is a generalization 
of the SK mean field spin glass to a model where inter- 
action is given by products of p spins (the p — 2 model 
coincides with SK): it turns out to be a promising candi- 
date to understand the physics of structural glasses 0,01- 
Numerical simulations of this model are extremely CPU 
time and memory consuming: one needs to store order 
of V p quenched random couplings, and the time needed 
to run a Metropolis sweep increases with volume like V p , 
as compared to V for short range models. Furthermore, 
glassy slow dynamics has to be studied for very large 
times. Here why study the p = 3 fully connected model: 
it is already quite expensive, but the cheapest of its class 
(we do not loose in generality since all models with p > 3 
are believed to have the same universal behavior). 

In the next section we give details about our numeri- 
cal simulations. We use parallel tempering Q, an opti- 
mized Monte Carlo method very effective on spin glass 
systems As an added bonus we present in the fol- 
lowing section very precise results about the statics of 
the system: these are by far the most accurate results 
obtained for equilibrium expectations in this model (see 
H for the former state of the art). They allow to com- 
pare the N = oc analytic results with finite N accurate 
values, and to get a feeling about what a "large" lattice 
really is. They also help in making us confident that we 
are really reaching thermodynamic equilibrium. The ev- 
idence collected here is interesting since, although this 
model is very important for reaching an understanding 
of the glassy phase, it has never been studied numerically 
with high accuracy, because of the extreme difficulty of 
the simulation. 

In a fourth section, we present our main results, about 
time scales in the model. We find that the Ioffe-Lopatin- 
Sherrington approach leads indeed to the correct esti- 
mate, and that ln(r Q ) oc V. We end by drawing our 
conclusions. 
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NUMERICAL SIMULATIONS 

The infinite range, fully connected, 3-spin spin glass is 
defined by the Hamiltonian 



H — — ^ JiJ,k <?i <7j &k , 
i<j<k 



(1) 



where the quenched random exchange coupling Jij t k gov- 
ern the interaction among triples of spins and can take 
one of the two values ±y3/V with probability one half. 
The choice we made of binary couplings allows far supe- 
rior performances of the computer code. 

The model has a complex phase diagram^, fiol Hl| . 
From a thermodynamic point of view it has three phases. 
At high T it is in a paramagnetic replica symmetric (RS) 
phase. At T = T c the system enters a one-step replica 
symmetry breaking (1RSB) phase, while at a lower tem- 
perature value Tq < T < T c it enters a third phase that 
will not concern us here. 

The order parameter is the usual overlap q. Its prob- 
ability distribution function can be written as 



(2) 



where, as usual, the brackets denote a thermal average 
and the over-line denotes an average over the quenched 
random couplings. In the RS phase, in the infinite vol- 
ume limit, two different replica have zero overlap. In the 
1RSB phase two different replica have overlap qo = 
with probability to, and overlap q\ > with probability 
1 — to (both qi and to depend on temperature, and should 
be written q\{T) and m{T)). The values of to and q\ are 
solutions of a set of two coupled integral equations. The 
condition to = 1 gives the value of the critical tempera- 
ture T c w 0.6514 (the expansion around the large p limit 
of gives T c « 0.6671, and T G » 0.24). It turns out 
that lim T ^ T - qi(T) ^ 0, namely the model has discon- 
tinuous 1RSB. 

Let us start by giving some details about our simula- 
tion. We study systems with V = 32, 64, 96, 128, 160 
and 192 spins (the needed CPU time increases with vol- 
ume like V 3 ). We first thermalize the system using the 
parallel tempering optimized Monte Carlo procedure 0, 
with a set of 15 T values in the range 0.3 — 1.2 (i.e. 
AT = 0.06) for the three smallest volumes and with 
AT = 0.03 in the temperature range 0.3 — 0.9 for the 
three largest volumes. We take advantage of the binary 
distribution of the couplings to use the multi-spin coding 
technique [12j, with a one order of magnitude gain in up- 
date speed. We perform 4 • 10 5 iterations (one iteration 
consists of one Metropolis sweep of all spins plus one tem- 
pering update cycle), and store the final well equilibrated 
configurations. 

We then, for studying the dynamical behavior of the 
system, start updating these equilibrium configurations 
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FIG. 1: The internal energy as a function of the temperature. 
The continuous curve is from the analytical solution of the 
model. Numerical data are for system sizes from 32 to 192. 
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FIG. 2: As in figure Q but for the specific heat. 
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FIG. 3: The order parameter q as a function of the tempera- 
ture. The continuous curve is from the analytical solution of 
the model. Data are for system sizes from 32 to 192. 
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FIG. 4: Probability distribution function of the overlap, for 
system sizes from 32 to 192. Here T = 0.42. 

with a simple Metropolis dynamics, and perform 2 • 10 6 
Metropolis sweeps. The number of disorder samples is 
1000 for the V = 32 system, 500 for V = 64 and 200 for 
the other, larger, sizes. As usual, the program simulates 
the independent evolution of two replicas, in order to 
compute the overlap q. All statistical error estimates are 
done using a jackknife analysis of the fluctuations among 
disorder samples, with 20 jackknife bins. 

We have used the second half of the 4 • 10 5 thermal 
sweeps, i.e. the last 2 • 10 5 , to measure static quantities. 
This is on the one side an interesting result on his own 
(since we can compare to the analytic, N = oo, result, 
and get hints about how finite N corrections work), and 
allows, on the other side, to check thcrmalization. No- 
tice that our largest thermalized system is more than five 
times larger than the largest analyzed before y| : we will 
present next these results, before discussing the dynam- 
ical behavior of the system. 
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FIG. 5: Binder parameter as a function of the temperature. 
The continuous curve is from the analytical solution of the 
model. 
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FIG. 6: The A parameter as a function of the temperature. 
The continuous curve is from the analytical solution of the 
model. 



STATIC BEHAVIOR 

We show in figure ^ our data for the internal energy 
E = (H)/V as a function of the temperature (we av- 
erage this single replica quantity over the two copies of 
the system that we follow in parallel), together with the 
analytical result 

E=-±(l-(l-m)ql) . (3) 

The finite volume numerical data converge nicely toward 
the infinite volume limit analytical curve. This statement 
can be make quantitative by fitting to some assumed an- 
alytical finite size behavior. Taking T = 0.3 as an exam- 
ple, a good representation of the data is obtained with 
the simple form Ey = Eoo + aV^ 1 /' : by using for Eoo 
the exact analytic value, one finds that a best fit gives 
i9 = 1.13 ± 0.01. In figure [3 we show the specific heat: 



again, the approach to the infinite volume result is very 
clear in the numerical data. 

We show in figure |3| the order parameter q = (q), to- 
gether with the theoretical result (1 — m)qi, as a function 
of the temperature: we find again an excellent conver- 
gence to the infinite volume limit. 

We show in figure^the overlap probability distribution 
(for T — 0.42). Data are in agreement with the predic- 
tion of a bimodal probability distribution, that becomes 
sharper and sharper as V grows: one peak is around zero 
and the other around some positive value of q. Between 
the two peaks the probability distribution vanishes in the 
large volume limit, as predicted by the 1RSB picture (and 
in contrast to ooRSB). 

Figures and El are for the usual Binder parameter 

*4(3-Jg), (4) 
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and for the A parameter that measures order-parameter 
fluctuations 



.4 



(I 2 )' 



(5) 



As shown in reference 8], B and A have very simple 
expression as a function of to, since for q = there is 
no residual dependence on qi and one gets that both B 
and A are equal to zero in the high T phase, and that: 



B = 
A = 



2 -3m 

2(1 - to) 
m 

3(1 -to) 



T G <T < T c , 
T G <T<T C . 



(G) 
(7) 



Note that we are in a non-standard case where some di- 
mensionless quantities diverge at T c : in this model both 
B and A diverge as T — > T~ . This is in marked contrast 
with usual cases, where both B and A have a finite limit 
for all temperatures as V — > oo. 

Our numerical data are consistent with the predictions: 
namely they show a zero limit in the RS phase, and a 
nonzero limiting curve, that diverges as T — > T c , in the 
1RSB phase (data are consistent with both the maximum 
of A and the minimum of B being proportional to V) . In 
other words, in this situation one finds, neither for B nor 
for A, a fixed point where the curves for increasing values 
of V cross. Absence of crossing for the Binder parameter 
has been also observed in for the infinite range 3 state 
Potts spin glass model (a model with continuous 1RSB) 
and in |j| for the 10 state model (a model with discon- 
tinuous 1RSB). The fact that A has a non-trivial limit 
as V grows shows clearly that the model has non-zero 
order parameter fluctuations (OPF) in the 1RSB phase, 
and the effectiveness of the parameter A to determine 
whether OPF holds or not. 



DYNAMICAL BEHAVIOR 

The main result of this note is the precise quantitative 
measurement of the typical time scales of the model and 
of their scaling behavior. As we have already discussed 
we use the approach of Q, Q , considering the equilibrium 
dynamics of thermalized configurations. 

We measure the time dependent overlap of two spin 
configurations: 



9(0, t) 



-^5>i(0>,(t) 



(8) 



where we start, at time zero, from an equilibrium spin 
configuration (see the former two sections), and the time 
t is measured in units of Metropolis sweeps (in this 
phase of the numerical simulation we are using the sim- 
ple Metropolis dynamics: we use parallel tempering only 
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FIG. 7: q c (0,t) as a function of t (in log scale) in units of 
Metropolis sweeps. Here T = 0.48. 
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FIG. 8: tx (in log scale) as a function of the number of spins 
V, for X = 0.5, 0.55, 0.6 and 0.65. T = 0.48. 



in the first phase of thermalization and equilibrium anal- 
ysis). The thermodynamic average (. . .) is taken by av- 
eraging over several initial configurations. As noticed in 
P| , the best and most effective way to proceed is to use a 
new disorder sample for every simulation: we do never re- 
peat simulations with different starting points in a given 
disorder sample, but always use new quenched random 
couplings when starting a new numerical run. This pro- 
cedure correctly averages over the thermal noise and over 
the random couplings, by minimizing the statistical in- 
certitude (mainly connected to the sample average): no 
bias is introduced. In what follows, the average (. . .) 
is accordingly an average over two independent replicas 
only. 

As t goes to infinity, q(Q, t) coincides with the static 
overlap, and it is distributed according to the probability 
distribution P(q). Accordingly (on a finite system) 



q c (0,t) = q(0,t) - q -> as t -► oo. (9) 
Our results for q c (0, t) as function of the Metropolis time 
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FIG. 9: Ti/2 (in log scale) as a function of the number of spins 
V. T — 0.54. 



(note the logarithmic scale) for T = 0.48 can be found 
in figured [T c ~ 0.6514, so we are at T ~ §T C ). The 
behavior of <? c (0,i) is very clear. On the largest lattices 
we are able to study, the time evolution has two succes- 
sive regimes as t grows: first g c (0,i) decays very slowly 
toward the value q± (the equivalent of qEA in spin glass 
models). It stays in this first regime longer and longer as 
the system size increases (and would stay there forever in 
the infinite system size limit). In a second regime, q c (0, t) 
goes to zero, and the system equilibrates. 

We need now to define a typical time scale of the equi- 
librium dynamics. We use the simple approach based on 
selecting a time such that the (normalized) correlation 
decreases below a given threshold X. We tried several 
values of this threshold to check the stability of this def- 
inition. So our time scale Tx is defined as the time such 
that q c (Q, t — tx) — X q c (0, 0). Technically if the numer- 
ical data cross the level X more than once we average over 
all crossing point (this is not a big issue, since we never 
find a real ambiguity, but only sometimes a small wig- 
gling, due to statistical fluctuations, very local in time). 
We now examine the values of Inrx as a function of V. 
If there is only one diverging time scale in the problem, 
the behavior of tx as function of V should not depend on 
X. We show the data obtained by this analysis in figure 
|H1 (for T — 0.48): there is indeed a remarkable level of 
universality, and the scaling law of the data is clearly in- 
dependent of X, making us confident about the quality of 
our approach. Notice that we are changing the threshold 
X in a large range of values, and the results are stable: 
this is a very good indication toward the fact that we are 
really determining the physical time scales. 

In figure 01 we show our data for t\i% at T = 0.54, 
together with the best fit of the data to the expression 

r{V)=Ae BV \ (10) 

with the values A = 35 ± 10, B = (4.5 ± 2.0) ■ 10~ 2 



e = 1.0 ± 0.1, with a value of chi-square \ 2 — 1-8 per 
degree of freedom. 

These numerical results are in remarkable agreement 
with the predictions of Q : in the broken symmetry phase 
the typical and relevant time scales grow here according 
to ln(r) oc V, as opposed to the Sherrington-Kirkpatrick, 
p = 2 model, where they grow according to ln(r) oc 
[2j. Here the absence of visible sub-leading corrections 
is remarkable, in contrast again to the situation of the 
Sherrington-Kirkpatrick model where the small expo- 
nent 1/3 makes things difficult. 



CONCLUSIONS 

We have studied the 3-spin infinite range spin glass 
model. We have determined with good accuracy equilib- 
rium properties on lattices of reasonable size, gaining in 
this way an accurate control of finite size effects. The 
main point of this note has been to investigate the equi- 
librium dynamics of the model, and to establish that the 
time scales of the model grow with system size V ac- 
cording to the law lnr a oc V: this is in nice agreement 
with the prediction of the theoretical approach by 0. 
Technically, from the computational point of view, these 
results are worthy because simulating a 3-spin infinite 
range model is a very non-trivial task: here the com- 
puter time needed to perform a full update of the lattice 
increases as V 3 for large volumes. We have succeeded on 
the one side to study the statics of the model on lattice 
more than five times bigger than the largest ones used 
before, and on the other side to determine with good ac- 
curacy severely increasing time scales: we consider that 
as a worthy achievement. 
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